function irfs = irf_ext_instr(bet_ident, Y, T, S, p)

% Compute IRFs for a monetary policy shock
% betas is the matrix of estimated coefficients, Y is the dependent vector
% of observations and T is the length of the irf. S is the identified
% structural shock response

shock = [1; zeros(size(Y,2)-1,1)]; %shock to the interest rate

% Eliminate constants and exovars from the original VAR coefficients
betas = bet_ident(1:end-1,:);

irfs  = zeros(size(Y,2),T+p);

% Simulate from p+1 onwards
for tt = 1+p:T+p
    if tt == p+1
        irfs(:,tt) = [S,zeros(size(Y,2),size(Y,2)-1)]*shock;
    else
        xx          = fliplr(irfs(:,tt-p:tt-1));
        irfs(:,tt)  = betas'*reshape(xx, size(xx,1)*size(xx,2),1); 
    end
end

%Eliminate the zeros before p
irfs = irfs(:,p+1:end);

end
